library(dplyr)
library(tidyr)
library(ggplot2)
library(rgdal)
library(raster)
setwd("~/Desktop/r_geospatial_102219/")
Error in setwd("~/Desktop/r_geospatial_102219/") :
cannot change working directory
GDALinfo("./data/HARV_dsmCrop.tif")
rows 1367
columns 1697
bands 1
lower left origin.x 731453
lower left origin.y 4712471
res.x 1
res.y 1
ysign -1
oblique.x 0
oblique.y 0
driver GTiff
projection +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
file ./data/HARV_dsmCrop.tif
apparent band summary:
GDType hasNoDataValue NoDataValue blockSize1 blockSize2
1 Float64 TRUE -9999 1 1697
apparent band statistics:
Bmin Bmax Bmean Bsd
1 305.07 416.07 359.8531 17.83169
Metadata:
AREA_OR_POINT=Area
Note: raster() will create a single-band raster. stack() and brick() can be used to load/create a multi-band raster (more on that later)
harv_dsm <- raster("./data/HARV_dsmCrop.tif")
harv_dsm
class : RasterLayer
dimensions : 1367, 1697, 2319799 (nrow, ncol, ncell)
resolution : 1, 1 (x, y)
extent : 731453, 733150, 4712471, 4713838 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : /Users/francesdavenport/Documents/teaching/r_geospatial/r-geospatial-102219/data/HARV_dsmCrop.tif
names : HARV_dsmCrop
values : 305.07, 416.07 (min, max)
plot(harv_dsm)
harv_dsm_df <- as.data.frame(harv_dsm, xy = TRUE)
str(harv_dsm_df)
'data.frame': 2319799 obs. of 3 variables:
$ x : num 731454 731454 731456 731456 731458 ...
$ y : num 4713838 4713838 4713838 4713838 4713838 ...
$ HARV_dsmCrop: num 409 408 407 407 409 ...
ggplot() +
geom_raster(data = harv_dsm_df, aes(x, y, fill = HARV_dsmCrop)) +
scale_fill_viridis_c()
harv_dtm <- raster("./data/HARV_dtmCrop.tif")
plot(harv_dtm)
Let’s calculate canopy height by subtracting the DTM from the DSM
harv_canopy <- harv_dsm - harv_dtm
plot(harv_canopy)
harv_canopy <- overlay(harv_dsm, harv_dtm,
fun = function(r1, r2) { return( r1 - r2) })
plot(harv_canopy)
Let’s calculate the mean height in our DTM and subtract that from our DTM
harv_dtm2 <- harv_dtm - cellStats(harv_dtm, stat='mean')
plot(harv_dtm2)
harv_dtm3 <- calc(harv_dtm, function(x) ifelse(x > 340, x, NA))
plot(harv_dtm3)
Let’s write a GeoTiff file. R will recognize the desired output file type based on the file extension (see documentation for options)
writeRaster(harv_canopy, "./data/HARV_canopy_height.tif")
Using the sp package:
aoi_sp <- shapefile("./data/HarClip_UTMZ18.shp")
aoi_sp
class : SpatialPolygonsDataFrame
features : 1
extent : 732128, 732251.1, 4713209, 4713359 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
variables : 1
names : id
value : 1
plot(aoi_sp)
library(sf)
Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
aoi <- st_read("./data/HarClip_UTMZ18.shp")
Reading layer `HarClip_UTMZ18' from data source `/Users/francesdavenport/Documents/teaching/r_geospatial/r-geospatial-102219/data/HarClip_UTMZ18.shp' using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 1 field
geometry type: POLYGON
dimension: XY
bbox: xmin: 732128 ymin: 4713209 xmax: 732251.1 ymax: 4713359
epsg (SRID): 32618
proj4string: +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
plot(aoi)
ne_states <- st_read("./data/Boundary-US-State-NEast.shp")
Reading layer `Boundary-US-State-NEast' from data source `/Users/francesdavenport/Documents/teaching/r_geospatial/r-geospatial-102219/data/Boundary-US-State-NEast.shp' using driver `ESRI Shapefile'
Simple feature collection with 12 features and 9 fields
geometry type: MULTIPOLYGON
dimension: XYZ
bbox: xmin: -80.51989 ymin: 37.91685 xmax: -66.9499 ymax: 47.45716
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
head(ne_states)
Simple feature collection with 6 features and 9 fields
geometry type: MULTIPOLYGON
dimension: XYZ
bbox: xmin: -79.76212 ymin: 37.91685 xmax: -66.9499 ymax: 47.45716
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
STATEFP STATENS AFFGEOID GEOID STUSPS NAME LSAD ALAND AWATER
1 11 01702382 0400000US11 11 DC District of Columbia 00 158350578 18633500
2 24 01714934 0400000US24 24 MD Maryland 00 25147575220 6983455225
3 36 01779796 0400000US36 36 NY New York 00 122054577774 19242052501
4 09 01779780 0400000US09 09 CT Connecticut 00 12542396439 1814978794
5 23 01779787 0400000US23 23 ME Maine 00 79885244456 11749091526
6 25 00606926 0400000US25 25 MA Massachusetts 00 20203936287 7131805598
geometry
1 MULTIPOLYGON Z (((-77.11976...
2 MULTIPOLYGON Z (((-76.04621...
3 MULTIPOLYGON Z (((-72.01893...
4 MULTIPOLYGON Z (((-73.69594...
5 MULTIPOLYGON Z (((-68.92401...
6 MULTIPOLYGON Z (((-70.27553...
plot(ne_states)
names(ne_states)
[1] "STATEFP" "STATENS" "AFFGEOID" "GEOID" "STUSPS" "NAME" "LSAD" "ALAND" "AWATER" "geometry"
ggplot() +
geom_sf(data = ne_states, aes(fill = NAME))
st_bbox(ne_states)
xmin ymin xmax ymax
-80.51989 37.91685 -66.94989 47.45716
st_crs(ne_states)
Coordinate Reference System:
EPSG: 4326
proj4string: "+proj=longlat +datum=WGS84 +no_defs"
nrow(ne_states)
[1] 12
mass <- ne_states %>%
filter(NAME == "Massachusetts")
ggplot() +
geom_sf(data = mass, aes(fill = NAME))
In addition to polygon shapefiles, we can read point and line data:
tower <- st_read("./data/HARVtower_UTM18N.shp")
Reading layer `HARVtower_UTM18N' from data source `/Users/francesdavenport/Documents/teaching/r_geospatial/r-geospatial-102219/data/HARVtower_UTM18N.shp' using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 14 fields
geometry type: POINT
dimension: XY
bbox: xmin: 732183.2 ymin: 4713265 xmax: 732183.2 ymax: 4713265
epsg (SRID): 32618
proj4string: +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
roads <- st_read("./data/HARV_roads.shp")
Reading layer `HARV_roads' from data source `/Users/francesdavenport/Documents/teaching/r_geospatial/r-geospatial-102219/data/HARV_roads.shp' using driver `ESRI Shapefile'
Simple feature collection with 13 features and 15 fields
geometry type: MULTILINESTRING
dimension: XY
bbox: xmin: 730741.2 ymin: 4711942 xmax: 733295.5 ymax: 4714260
epsg (SRID): 32618
proj4string: +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
site_map <- ggplot(roads) +
geom_sf() +
geom_sf(data = tower, col = "red", size = 10, shape = "*")
Note: ggplot is displaying coordinates as lat/lon, even though our vector data has the UTM projection!
st_crs(tower)
Coordinate Reference System:
EPSG: 32618
proj4string: "+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs"
tower_buffer <- st_buffer(tower, 200)
site_map +
geom_sf(data = tower_buffer, fill = "green", alpha = 0.2)
Note: We could use crop to clip the raster to the bounding box of an object, but not it’s actual boundary
canopy_masked <- mask(harv_canopy, tower_buffer)
canopy_masked
class : RasterLayer
dimensions : 1367, 1697, 2319799 (nrow, ncol, ncell)
resolution : 1, 1 (x, y)
extent : 731453, 733150, 4712471, 4713838 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : layer
values : 0, 38.16998 (min, max)
plot(canopy_masked)
buffer_canopy_values <- extract(harv_canopy, tower_buffer)
Extract returns a list of vectors (one for each polygon object), but we only have one object so we can unlist the result
buffer_canopy_values <- unlist(buffer_canopy_values)
hist(buffer_canopy_values)
ndvi_files <- list.files("./data/NDVI/",
full.names = TRUE,
pattern = ".tif$")
harv_ndvi <- stack(ndvi_files)
harv_ndvi
class : RasterStack
dimensions : 5, 4, 20, 13 (nrow, ncol, ncell, nlayers)
resolution : 30, 30 (x, y)
extent : 239415, 239535, 4714215, 4714365 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=19 +ellps=WGS84 +units=m +no_defs
names : X005_HARV_ndvi_crop, X037_HARV_ndvi_crop, X085_HARV_ndvi_crop, X133_HARV_ndvi_crop, X181_HARV_ndvi_crop, X197_HARV_ndvi_crop, X213_HARV_ndvi_crop, X229_HARV_ndvi_crop, X245_HARV_ndvi_crop, X261_HARV_ndvi_crop, X277_HARV_ndvi_crop, X293_HARV_ndvi_crop, X309_HARV_ndvi_crop
min values : 2732, 1534, 1917, 5669, 8685, 8768, 8633, 8686, 8297, 7491, 277, 484, 4396
max values : 5545, 3736, 3600, 6394, 8903, 9140, 8945, 8967, 8651, 8607, 366, 661, 6359
crs(harv_ndvi)
CRS arguments:
+proj=utm +zone=19 +ellps=WGS84 +units=m +no_defs
crs(harv_dtm)
CRS arguments:
+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
harv_ndvi <- projectRaster(harv_ndvi, crs = crs(harv_dtm))
no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf
ndvi_mean <- cellStats(harv_ndvi, mean)
plot(ndvi_mean)
We could do this by finding the cell overlapping with the tower and extracting the value at that point
tower_xy <- st_coordinates(tower)
tower_row <- rowFromY(harv_ndvi, tower_xy[2])
tower_col <- colFromX(harv_ndvi, tower_xy[1])
tower_ndvi <- getValues(harv_ndvi, row = tower_row)[tower_col,]
plot(tower_ndvi)
tower_ndvi <- raster::extract(harv_ndvi, tower)
buffer_ndvi <- raster::extract(harv_ndvi, tower_buffer, fun = mean)
buffer_ndvi
X005_HARV_ndvi_crop X037_HARV_ndvi_crop X085_HARV_ndvi_crop X133_HARV_ndvi_crop X181_HARV_ndvi_crop
[1,] 3709.251 2507.582 2571.568 6004.585 8780.76
X197_HARV_ndvi_crop X213_HARV_ndvi_crop X229_HARV_ndvi_crop X245_HARV_ndvi_crop X261_HARV_ndvi_crop
[1,] 8941.279 8783.959 8818.749 8503.574 7955.772
X277_HARV_ndvi_crop X293_HARV_ndvi_crop X309_HARV_ndvi_crop
[1,] 332.0345 561.2631 5457.627
ndvi_ts <- as.data.frame(unlist(buffer_ndvi)[1,])
names(ndvi_ts) <- "meanNDVI"
julian_day <- gsub("X|_HARV_ndvi_crop", "", row.names(ndvi_ts))
origin <- as.Date("2011-01-01")
ndvi_ts$date <- origin + as.integer(julian_day) - 1
ndvi_ts
meanNDVI date
X005_HARV_ndvi_crop 3709.2514 2011-01-05
X037_HARV_ndvi_crop 2507.5816 2011-02-06
X085_HARV_ndvi_crop 2571.5682 2011-03-26
X133_HARV_ndvi_crop 6004.5853 2011-05-13
X181_HARV_ndvi_crop 8780.7595 2011-06-30
X197_HARV_ndvi_crop 8941.2789 2011-07-16
X213_HARV_ndvi_crop 8783.9587 2011-08-01
X229_HARV_ndvi_crop 8818.7488 2011-08-17
X245_HARV_ndvi_crop 8503.5736 2011-09-02
X261_HARV_ndvi_crop 7955.7723 2011-09-18
X277_HARV_ndvi_crop 332.0345 2011-10-04
X293_HARV_ndvi_crop 561.2631 2011-10-20
X309_HARV_ndvi_crop 5457.6267 2011-11-05
ggplot(ndvi_ts, aes(date, meanNDVI)) +
geom_point() +
geom_line()
Remember, our vector data includes roads and our area of interest, in addition to making our own areas using things like st_buffer! Try to plot your results and discuss it with your neighbor